
#install.packages("mediation")
library("mediation")

set.seed(12345)

dat = read.csv("Mediation.csv", header=T)

for (i in 1:ncol(dat)) {
  dat[,i] = as.numeric(dat[,i])
}


mod.med.bonding = lm(disagreement ~ partybonding + 
             female + age + education + income + rel2 + rel3 + rel4 + nonwhite + pidstrength + repdummy, data=dat)

mod.med.bridging = lm(disagreement ~  partybridging +
               female + age + education + income + rel2 + rel3 + rel4 + nonwhite + pidstrength + repdummy, data=dat)

mod.dv.interaction = lm(interaction ~ disagreement + partybonding + partybridging +
                        female + age + education + income + rel2 + rel3 + rel4 + nonwhite + pidstrength + repdummy, data=dat)

mod.dv.interest = lm(polinterest ~ disagreement + partybonding + partybridging +
                     female + age + education + income + rel2 + rel3 + rel4 + nonwhite + pidstrength + repdummy, data=dat)

mod.dv.talk = lm(politicaltalk ~ disagreement + partybonding + partybridging +
                 female + age + education + income + rel2 + rel3 + rel4 + nonwhite + pidstrength + repdummy, data=dat)



###########################################
## Calculation of Mediation Effect
###########################################

med.interaction.bonding = mediate(mod.med.bonding, mod.dv.interaction, treat="partybonding", mediator="disagreement")
med.interaction.bridging = mediate(mod.med.bridging, mod.dv.interaction, treat="partybridging", mediator="disagreement")

med.interest.bonding = mediate(mod.med.bonding, mod.dv.interest, treat="partybonding", mediator="disagreement")
med.interest.bridging = mediate(mod.med.bridging, mod.dv.interest, treat="partybridging", mediator="disagreement")

med.talk.bonding = mediate(mod.med.bonding, mod.dv.talk, treat="partybonding", mediator="disagreement")
med.talk.bridging = mediate(mod.med.bridging, mod.dv.talk, treat="partybridging", mediator="disagreement")


sink("Mediation ACME.txt")
print("============ Interaction Bonding ============")
summary(med.interaction.bonding)
print("============ Interaction Bridging ============")
summary(med.interaction.bridging)

print("============ Interest Bonding ============")
summary(med.interest.bonding)
print("============ Interest Bridging ============")
summary(med.interest.bridging)

print("============ Talk Bonding ============")
summary(med.talk.bonding)
print("============ Talk Bridging ============")
summary(med.talk.bridging)
sink(NULL)


###########################################
## Sensitivity Analysis
###########################################

width = 3200
height = 2400
res = 300

medsens.interaction.bonding = medsens(med.interaction.bonding)
medsens.interaction.bridging = medsens(med.interaction.bridging)

medsens.talk.bonding = medsens(med.talk.bonding)
medsens.talk.bridging = medsens(med.talk.bridging)


png(filename = "Sensitivity Interaction Bonding.png", width=width, height=height, res=res)
plot.medsens(medsens.interaction.bonding, main="Outparty Interaction - Bonding")   
axis(side=1, at=seq(-1,1,by=.1))
dev.off()

png(filename = "Sensitivity Interaction Bridging.png", width=width, height=height, res=res)
plot.medsens(medsens.interaction.bridging, main="Outparty Interaction - Bridging")   
axis(side=1, at=seq(-1,1,by=.1))
dev.off()

png(filename = "Sensitivity Talk Bonding.png", width=width, height=height, res=res)
plot.medsens(medsens.talk.bonding, main="Political Discussion - Bonding")   
axis(side=1, at=seq(-1,1,by=.1))
dev.off()

png(filename = "Sensitivity Talk Bridging.png", width=width, height=height, res=res)
plot.medsens(medsens.talk.bridging, main="Political Discussion - Bridging")   
axis(side=1, at=seq(-1,1,by=.1))
dev.off()


interaction = data.frame(interaction.bo.lower = medsens.interaction.bonding$lower.d0,
                         interaction.bo.upper = medsens.interaction.bonding$upper.d0,
                         interaction.br.lower = medsens.interaction.bridging$lower.d0,
                         interaction.br.upper = medsens.interaction.bridging$upper.d0
                         )

talk = data.frame(talk.bo.upper = medsens.talk.bonding$upper.d0,
                  talk.br.lower = medsens.talk.bridging$lower.d0,
                  talk.br.upper = medsens.talk.bridging$upper.d0
                      )
out = cbind(seq(-.9,.9,by=.1), interaction, talk)
colnames(out)[1] = "rho"

write.csv(out, file="Mediation Sensitivity.csv", row.names=FALSE)
